Density of states determined from Monte Carlo simulations 
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We describe method for calculating the density of states by combining several canonical monte 
carlo runs. We discuss how critical properties reveal themselves in g(e) and demonstrate this by 
applying the method several different phase transitions. We also demonstrate how this can used 
to calculate the conformal charge, where the dominating numerical method has traditionally been 
transfer matrix. 
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I. INTRODUCTION 

Since it was devised by Metropolis in 19530 Monte 
Carlo (MC) simulations based on the Metropolis algo- 
rithm have had a tremendous impact on physics; SIAM 
recently rated the algorithm among the ten most influen- 
tial numerical algorithms of the previous century 0- For 
a historical summary, and a review of modern MC meth- 
ods we refer to the proceedings of the conference hosted 
to celebrate the 50th anniversary of the algorithm^. For 
general references to Monte Carlo simulations see e.g. 
Refs. SHI 

The Metropolis algorithm is well suited to calculate 
quantities which can be expressed as 
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i.e. as averages of values obtained by operating an op- 
erator O on a series of states \ipi). Focusing on phase 
space Eq.^can be denoted a local estimator, in the sense 
that only one point in phase space is involved at a time. 
Some quantities like entropy and free energy can not be 
expressed like Eq. ^ their evaluation requires simultane- 
ous knowledge of global portions of phase space. Entropy 
and free energy can in principle be obtained by thermo- 
dynamic integration^]]], 



F(T) = U{T) -T [ dT 
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but this technique does not seem to be much used. 

Eq. ^ represents the absolutely simplest way to get 
MC results. A simulation produces a series of states 
distributed according to the Boltzmann distribution, the 
mean over these states is calculated. Both the initial step 
of obtaining the data, and the final post-processingcan 
be done differently. With multicanonical sampling 0> S @ 
the Markov chain is altered to (ideally) yield a flat energy 
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histogram; and the results reweighted back afterwards. 
The Wang-Landau histogram method[l3,[l3 can be seen 
as combined data collection and post-processing; when 
the simulation is complete we have built up an estimate 
g(e) of the density of states (DOS). For some situations 
like first order transitions and disordered media these 
methods have very efficient. 

During the simulation we can build up an estimate of 
the complete density Pfs(0), and clearly it would be bene- 
ficial to utilize this information. This insight is the ke y to 
histogram methods. In 1989 Ferrenberg and Swendsen|l2j 
published a method to combine results obtained at dif- 
ferent couplings. The method was highly efficient, and 
Ferrenberg-Swendsen reweighting has become an essen- 
tial tool for MC practitioners. The use of rawdata from 
several couplings allow for reweighting to a much broader 
range of couplings than ordinary single histogram meth- 
ods. 

In 1990 Alves, Berg and Villanova (ABV)[I3 de- 
veloped a variation of the Ferrenberg-Swendsen multi- 
histogram technique; specifically targeted at calculating 
the density of states. To apply the FS method one must 
solve a set of nonlinear equations self consistently, this 
can fail if the overlap between the various histograms is 
insufficient. This is not the case for the ABV method 
which can always be applied as long as every histogram 
has finite overlap with at least one other histogram. We 
have developed a method to calculate DOS which is a 
minor variation of ABV's original method. 

The density of states is an elusive quantity, and not 
very much used in statistical mechanics. In addition 
to presenting a method to calculate g(e) we have there- 
for also briefly discussed statistical mechanics based on 
g(e) in section ITTT1 and several applications in IIVI Some 
of these applications are well known results from tradi- 
tional canonical thermodynamics, however there are also 
properties which are more easily learned based on micro- 
canonical thermodynamics. 

The main focus of this paper is to determine the den- 
sity of states from canonical Monte Carlo simulation. 
The density of states is the central quantity in micro- 
canonical thermodynamics; hence this naturally becomes 
an important formalism for further analysis of the DOS 
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based results. The study of microcanonical thermo- 
dynamics has seen increasing interest the latest years; 
see e.g . Ref. for a general introduction, and Refs. 
ITU lid" for some recent applications. 

The rest of the paper is organised as follows: In section 
ITTI we present the algorithm to calculate the density of 
states. Section IIIII is devoted to a short discussion of 
statistical mechanics based on g(e). In the final section 
IIVI we use the algorithm to study several different phase 
transitions. 



II. ESTIMATING g(e) 



sity of states in Eq.0]has an index j3 to indicate that the 
histogram was recorded at this coupling, but it does not 
have any intrinsic temperature dependence. In principle 
Eq. 0] can be used to estimate g(e) regardless of tem- 
perature, however practically only a small energy range 
around (E)(T) will be sampled with a sufficiently high 
frequency. 

Although Eq. 0] is useless as an immediate estimator 
for g(e), it provides a basis for combining results from 
different couplings to an estimator g(e) which can be ap- 
plied over the complete energy range. Given N different 
histograms hi(E) recorded at the couplings /?i > (3% > 
■ ■• > Pisr, we can combine them as 



When doing a MC simulation with the Metropolis al- 
gorithm the probability to be in a state i/j with energy e^> 
is proportional to 



(3) 



If we record a histogram of energies from a simulation 
at coupling j3; we get a histogram hp(e) which is propor- 
tional to g(e)e~ l3<L . Multiplying this histogram with e^ £ 
we get something which is proportional to g(e), i.e. 
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is an estimator for 5(e). In Eq. 0]e^ is a dimensionless 
constant of proportionality to be determined. The den- 
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hi(e) 
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to obtain an estimator which is usable over the complete 
e range [min e /ij(e),max e h(e)]. In Eq.EJtUjfe) is a weight 
function, which denotes the weight ascribed to histogram 
i in the estimation of g(e). The constants are deter- 
mined by joining the various histograms. 

The algorithm we have applied to determine £j is to 
set £1 to an arbitrary value, and then compute £j>i by 
minimising 



N—l N 
i—1 e 
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In §/3i( e )- ln ff/3j(c) 



As indicated in Eq. 0the central principle is to minimise 
the pairwise difference between all the g(e) estimates, 
where the estimates g$ i (e) are given according to Eq. 0] 
Minimising \ 2 with respect to £j gives N — l linear equa- 
tions which can be solved by e.g. LU decomposition. 
The algorithm described by ABV uses a different weight 
u>i(e) and the coefficients £j are determined from a re- 
cursive procedure; ^+1 is given by & and a function of 
the overlap between histogram ft,*(e) and hi+-y(e) (i and 
i + 1 are not necessarily ordered according to coupling, 
see Ref. Il3| for details) — . Apart from these differences 
this algorithm coincides with the one by ABV. 

When the coefficients & have been determined we have 
all the coefficients £j>i expressed in terms of £1. For dis- 
crete models with a finite ground state degeneracy go we 
can determine £1 by requiring g(eo) — go, or alternatively 
if the total number of states is known, this can be used 
to normalize g(e). In section Hvl we will consider both 
discrete models were the complete normalisation can be 



achieved, and continuous models were £1 must be left 
undetermined. 

Use of Eq. to determine 17(e) is in principle quite 
straightforward, but in practice it is important to be care- 
ful to avoid numeric underflow or overflow in intermediate 
steps, in particular the implementation must ensure that 
only lng(e) is needed in actual computations. 



III. STATISTICAL MECHANICS FROM g(e) 

Knowledge of g(e) is in principle equivalent to knowl- 
edge of the partition function Z(f3), hence all the proper- 
ties of a system are contained in g(e), however g{e) does 
not have a very prominent role in modern statistical me- 
chanics. We will therefor express some important results 
based on g{e) in this section, examples/applications are 
given in section IIVI The definition of temperature in 
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statistical mechanics is given by 



dlng(e) 
<9e ' 
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From this we find that the fundamental requirement 
CVCZ 1 ) > is equivalent to d e 2 In g(e) < 0. The limiting 
value d e hig(e) = C is the signature of a phase transi- 
tion. A finite e range with d e In g(e) = C means that the 
temperature is unchanged for this e range, i.e. an indi- 
cation of a first order transition; actually, as we shall see 
in section llV Bl this is slightly more complicated. When 
the width of the of linear part of lng(e) diminishes the 
first order transition is weakened; until <9 2 In 17(e) = in a 
isolated point only, this is the manifestation of a critical 
point. If we differentiate Eq. 0with respect to T we find 
the function Cy(e) 



-(d e \ng(e)f 
de l n 5(e) 



(8) 



From Eq. [S] we see that the critical properties, and in 
particular the critical exponent a, must be related to how 
3 2 In 17(e) approaches zero. To infer a directly from the 
behaviour of g(e) close to e c is difficult, but if we make 
the size dependence of 17(e) explicit we can use finite size 
scalingQ, [l|l|2(j. At the (pseudo)critical point in a finite 
system, Cy scales as i, d + a / 1/ ^ The factor (<9 e In g(e)) 2 in 
Eq. is just equal to /3 2 , hence the critical properties 
must come from the second derivative 



\d*hig(e,L)\L d oc L 



— a/ V 
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In general d e \ng(e, L) will also have finite size effects, 
however for this only the deviation from the thermody- 
namic value will show critical scaling. 

For microcanonical systems the externally specified 
variable is e, and not T, and critical sca ling is governed 
by the difference e — e c |, see e.g. Ref. |2(j. Using this 
Hiiller and Pleimling have calculated the order parame- 
ter exponent (3 from microcanonical data from the two 
and three dimensional Ising mo del [ill . 

When we have g(e) we can easily calculate F(T) and 
P(e,T) 



F(T) = -Tln^ 5 (e)e-^, 

e 

ff (e)e-* 



P(e,T) 



/3e' 
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From P(e, T) we can easily calculate the internal energy, 
and all moments thereof. If we in addition to e sam- 
ple other operators like the magnetisation, we can use 
P(e, T) to find thermal averages of arbitrary operators, 



(12) 



IV. SOME APPLICATIONS 

In this section we present various applications of the 
method presented in the preceding sections. In section 
IIV Al the results are benchmarked against the 2D Ising 
model, where F( T) has been determined exactly even 
for finite systems \2)\. In section TlV Bl we investigate the 
way phase transitions reveal themselves in g(e). In sec- 
tion ^^21 we calculate the finite size corrections to the 
free energy in a cylindrical geometry. For conformally 
invariant systems this is universal 23:, |24| . and can 
be used to calculate the conformal charge. We have de- 
termined the conformal charge for the 2D Ising model, 
and the 2D Q = 3 Potts model. In section llV Dl we dis- 
cuss the problems related to models with a continuous 
energy distribution; and show that method is useful also 
for these systems, although less so. Finally in section 
IIV El we have applied the method to a large dataset ob- 
tained from a previous study of the full Ginzburg Landau 
(GL) model. 



A. Comparisons with exact results 

Due to Onsagers exact solution|25| the 2D Ising model 
has been one of the most used benchmarks in statistical 
physics. For a rectangular lattice with periodic boundary 
the model has been solved in closed form even for finite 
systems j2^|, this constitutes a very convenient bench- 
mark for our approach. We have performed simulations 
on a 32 x 32 system with periodic boundary conditions, 
and verified that within statistical error both F(T) and 
Cy(T) agree with the exact values, see figures El 
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FIG. 1: This figure shows the estimated value of F(T) as 
symbols, and the exact value from jgjf as a solid line. The 
small ticks on T axis indicate T values were simulations have 
been performed. In the inset the dashed curve shows the 
relative error of the estimated values, and the solid lines are 
± an estimated statistical error. 



In Eq. E| (0) t is the mean of O for a given value of e. 



The statistical errors errors indicated on figuresnand[21 
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FIG. 2: This figure shows essentially the same as Fig. Q but 
for the specific heat Cv{T). 

are calculated by performing ten completely independent 
simulations. 



B. The signature of phase transitions in g(e) 

As discussed in section ITTTl all critical properties must 
be present in g(e). In this section we will discuss the 
critical properties of the Q — 3 and Q = 10 Potts model. 
The first model has a continuous phase transition with 
a = 1/3 and v — 5/6 j2|j, i.e. ajv — 0.4, the second 
model has a first order transition. 

First we consider the Q = 3 model, for this model the 
goal is to determine the ratio a/u from g(e, L). According 
to Eq. lathis can be done by considering how In g(e, L) 
vanishes when approaching the critical energy e c . Fig- El 
shows L d d^\ng(e, L) for different system sizes, and we 
can see that peak approaches zero with increasing system 
size. 
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FIG. 3: L d dt In g(e, L) for different system sizes, in the limit 
L — > oo this vanishes as L~ a l v . 



In Fig. 0] we have plotted min \L d &l \ng(e, L)\. i.e. the 
magnitude of the peak value for the curves in Fig. |3 as 
a function of L. 
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FIG. 4: This shows finite size scaling of L d min \dl g(e, L)\ in 
a log-log plot, the dashed line has slope —ajv ~ —0.29. 

The dashed line in Fig.0|has slope of —ajv w —0.29; 
this is a significant deviation from the exact value ajv — 
0.40, however we feel that these results are sufficient to 
demonstrate that the critical properties, and in particu- 
lar the exponents a and v are contained in <?(e). There 
is clearly significant finite size effects in d e In g(e, L) also; 
including the factor d t h\g{e) gives the "improved" value 
ajv ps 0.35, however this can not contribute in the 
L — > oo limit and we have therefor not included this 
factor in Fig. Q] Finally Fig. 0] is based on the second 
derivative of a sampled quantity; hence it will clearly 
be difficult to determine with high precision. In conclu- 
sion it is definitely possible to infer the ratio ajv from 
the properties of g(e,L), but it is certainly not the most 
suitable way for high precision measurements. Finally we 
mention that the remaining critical exponents can not be 
obtained from g(e), their value is based on the explicit 
choice of fields to represent the critical state. 

Altough all thermodynamic information about a phase 
transition is contained in F(T), it is only for a first order 
transition, where DtF(T) is discontinous at T c , that the 
phase transition stands out in F(T). Fig. [S] shows F{T) 
for the strongly transition in the two dimensional Q = 10 
Potts model; a discontinuity in dxF(T) at T « 0.71 is 
easily spotted. 

For second order transitions we had to revert to FSS 
to infer critical properties from 17(e); in the case of first 
order transitions we can make quite powerful statements 
from g(e) alone. Given a first order transition between 
the pure states t\ and £2 the probability P(e,T) is bi- 
modal, with distinct peaks at the pure energy levels ei 
and €2- The mixed states with energy e\ < e < £2 are 
exponentially suppressed, to reproduce this behaviour we 
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must have 

ln ff ( ei + Ae)<ln 5 ( ei )+/31n^^|-^^. (13) 

for < Ae < 62 — e\. Hence In 17(e) must increase weaker 
than linearly in the immediate vicinity of ei and then sub- 
sequently stronger than linearly afterwards such that the 
relation j3 c = d e lng(ei) = d e g(e2) is satisfied. If we insist 
on d^lng(e) < also in the interval e\ < e < £2 In g(e) 
must have a cusp in this interval; however these energy 
levels correspond to states which are manifest not equi- 
librium so it might be too strict to require &j In g(e) < 
in this particular interval. An extensive discussion of the 
region e\ < e < 62 can be found in discussed in Ref. I27I 
The various details of 17(e) around a first order transition 
are illustrated in Fig. |SJ 
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FIG. 5: The free energy F(T) for the Q = 10 Potts model 
which has a strong first order transition. Although there is 
inevitably some finite size rounding, the cusp in this figure 
is quite clear. Fig. Q shows a similar figure for a continuous 
transition, this clearly smooth in comparison. 



C. Conformal charge 

It is well known that critical systems are scale invari- 
ant; in addition the critical systems often possess further 
symmetries like translational and rotational invariance. 
Together these operations form a group Q. Under quite 
mild restrictions, in particular finite length interactions, 
the system is actually invariant also under transforma- 
tions which vary in space, this means that Q is the con- 
formal group. In particular for d = 2 this is a very pow- 
erful result, and the application of Conformal Field The- 
ory (CFT) has lead to many exact results for the critical 
state[28|. 

Consider an infinitely long strip of width W, due to 
the finite width there will be finite size corrections to the 
free energy density. One of the most fundamental results 
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FIG. 6: Results from the Q = 10 Potts model. The up- 
permost panel shows ln<j(e), the two dashed lines indicate 
the two energy levels ei and £2- The central panel shows 
lngi(e) — [/3(e — ei) + lng(ei)], which clearly shows that In g(e) 
has small but significant deviations from perfect linearity in 
the range ei < e < £2- The bottom panel shows P(e, T c ), i.e. 
Eq. 1111 at the critical point, and we can clearly see how the 
depression in P(e,T c ) originates from the features in In 17(e). 
This figure very closely resembles Fig. 1 of Ref. l2?l 



from conformal invariance is that the leading finite size 
correction for this system is universalis 01 

In Eq. ^] fw is the free energy density of the strip, 
Jb is the bulk free energy density of an infinite system 
and c is the conformal charge or anomaly. The conformal 
charge is a dimensionless number which uniquely charac- 
terises a given universality class. In two dimensions the 
Ising model has c = 0.5 and the Q — 3 Potts model has 
c = 0.8H2. 

Eq. 1141 is a finite size scaling expression which should 
be very useful for numeric evaluation of c. However, since 
the use of Eq. ^] requires knowledge of the free energy 
MC has not been extensively used; see however for in- 
stance 29] for a numerical test of another Conformal 
Field Theory conjecture by MC methods. The numerical 
evaluation of c has been dominated by transfer matrix 
methods jS, see however Ref. |^ for a study of the 
Q = 3 Potts model very similar to the present one. 

Using the method presented in this paper we have 
calculated c for the Ising model and the Q — 3 Potts 
model. We have considered cylindrical systems of length 
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L and circumference W, we have considered W = 
{4, 5, 6, 8, 10, 12, 16}, and for each W we have used L = 
{W, 2W, AW, 8W, (WW)}, L = 16W was only considered 
for the Ising model. From this we have extrapolated to 
find 
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FIG. 7: The extrapolated limit fw from Eg. I15l for the Ising 
and Q = 3 Potts model. The conformal charge is given by the 
slope in the W~ 2 — > limit. The solid linxe is a least squares 
fit to Eq. 1141 and the dashed line comes from including an 
additional term aW~ 4 in the fit. For the Ising model the 
error bars are smaller than the symbol size. 

Plots of fw are shown in Fig. [7| and c has been de- 
termined from a least squares fit to these curves. The 
curves in Fig. show that there are quite significant cor- 
rections to the W~ 2 term for small W, and the numer- 
ical results c = 0.49 ± 0.07 for the Ising model was ob- 
tained by excluding all W < 8. Including an additional 
fourth order term[30j, and including all the results we get 
c = 0.55±0.06. Most of the computational resources were 
spent on the Ising model, and the Q — 3 Potts results are 
lower quality, the solid line corresponds to c = 0.86±0.19, 
this was obtained by retaining only the four largest W 
values. 



D. Continuous systems 

The method we have presented can to some extent also 
be applied to systems with a continuous energy distribu- 
tion, however for these systems the full normalisation of 
g(e) is difficult. The method used to normalise 17(e) for 
discrete systems so far require that (i) the ground-state is 
sampled, and (ii) that the histograms have sufficient over- 
lap. For a system with a truly continuous energy distri- 
bution sampling of the ground-state requires T = 0, and 
this will will generate histograms without overlap. Even 
for models with a very small energy gap, like e.g. the 



Z q 32] model for large q a large fraction of the computa- 
tional resources must be spent close to the ground state 
to ensure that both the conditions are met. Attempts 
to generalise the Wang Landau histogram sampling to 
continuous systems are faced with essentially the same 
problems |33|. 

Due to the problems with normalisation we must gen- 
erally be content with a function In g(e) — 1115(e) + C 
where C is an unknown, dimensionless constant. This 
will induce a linear error AF(T) = — TlnC in the free 
energy, but since Eq. 1111 is independent of C all remain- 
ing thermodynamics will be unaffected by the incomplete 
normalisation. 

The Z q clock model is a planar spin model where the 
real angle <fi £ [0, 2tt) is approximated with the discrete 
variable Qi — i2n/q, in the limit q — > co the converges to 
the XY model. Numerical simulations of the XY model 
are customarily done using the Z q model with q "large 
enough" , values of q — 32 or 64 are often used. Further- 
more it has been shown that already at q = 5 the critical 
properties are governed by the XY critical point [32ll3^ |. 
To learn about the behaviour g(e) for continuous systems 
we have done a short simulation of a 32 x 32 system for 
the Z q model with q — 2048, according to the discussion 
above this should with a good capture the properties of 
the continuous q — > 00 system. Fig. [H] shows lng(e) for 
the Z q model with q = 2048 along with the Ising model. 
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FIG. 8: Comparison of In g(e) for the q — 2048 Z q model and 
the Ising model. Observe that the curve for the Z q model has 
been vertically shifted by an unknown amount, see text. 

For the Z q model the lowest lying of the sampled states 
has energy e_ , hence the Z q curve in Fig. |H1 terminates 
at this e value. This means that <?(e) is undetermined 
in the interval [eo,e~), where eo = —2L d is the ground- 
state energy. Furthermore overall normalisation is im- 
possible to determine, and we have just arbitrarily fixed 

Vg(e-) = 0. 

According to Eq. 1111 internal energy and specific heat 
depend only on the shape of lng(e), and not possible 
vertical offset. Combined with the knowledge from pre- 
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vious simulations: that the Z q model for large q cor- 
rectly captures the thermodynamics of the XY model, 
we can conclude that apart from the vertical offset Fig. El 
is a faithful representation of g(e) for the continuous XY 
model. Looking at Fig. |S|the most striking features are 
(i) g(e) is orders of magnitude larger for the Z q model 
than the Ising model, and (ii) the steep sloop of 1115(e) 
for the Z q model; actually any gap-less model must have 
lim e ^ eo d e g(e) — > 00 to produce a finite value for Cy(T) 
in the limit T — » 0. 

In section HVEl we will reanalyse a real dataset from 
a large scale simulations of the Ginzburg Landau model, 

I 




In Eq. ^3 <f> is a complex condensate field, A is the 
electromagnetic gauge field, and x and y are parame- 
ters. The system is driven through a phase transition 
by the temperature like parameter y, and the qualitative 
behaviour at the phase transition is governed by x. For 
large x amplitude fluctuations in 4> are suppressed, leav- 
ing only phase fluctuations, and the transition is second 
order. In the limit x — > the amplitude fluctuations 
dominate, and the transition is first order. For interme- 
diate values of x all the degrees of freedom influence the 
dynamics, and for x — xj the transition changes order 
at a tricritical point |35j. 

In 2001 we determined the tricritical value xj from 
large scale Monte Carlo simulations. Here we have re- 
analysed some of the data from this simulation. Fig. [5] 
shows results for the GL model for x < xj, i.e. first 
order transition. The upper panel of this figure can be 
compared with Fig. [S] and the two lower panels can be 
compared with the two upper panels of Fig. 

Fig. Ellis similar to Fig.^l but for x > xj, i.e. a second 
order transition. 
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this constitutes a real example of a continuous system. 



E. Reanalyzing Ginzburg Landau results 

The Ginzburg Landau (GL) model is one of the most 
studied models physics, and it is applied as 'meta-model' 
in a wide range of fields, see Ref. [35| for an extensive list 
of applications. In dimcnsionlcss form, the continuum 
version of the model is given by the functional integral 



(16) 

I 

Due to the difficulties mentioned in section llV Dl we are 
not able to calculate the overall normalisation go of g(e), 
nevertheless g(e) shows the critical behaviour discussed in 
section llV Bl and the behaviour of F(y) clearly separates 
between first and second order transitions. F 

In conclusion we feel that in the application to the GL 
model the method has proved itself, and furthermore that 
it provides interesting information even tough g$ can not 
be determined. 

Software to go through the steps described in sections 
Inland lTTTl can be down-loaded as a C library from the au- 
thors web-site http://www.ii.uib.no/~hove/libdos/ 
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